visualize_means <- function(.data,
                            .sample = c("all_respondents", "by_party", "did", "panel")){
  
  # 1. Depending on the sample visualization (either all, DiD, or panel),
  # set initial variables and legend scales for each visualization
  
  # All Respondents
  if(.sample == "all_respondents" | .sample == "by_party"){
    .color = sym("Party")
    scale_color = c( "#0015BC", "#DE0100")

    .x = "Estimated outcome"
    out0 <- .data
    
    position_scales <- list(scale_x_continuous(breaks = seq(0, 40, 10), limits = c(0, 40)),
                            scale_x_continuous(breaks = seq(0.6, 1, 0.1), limits = c(0.6, 1)))
  # DiD Analysis
  } else if (.sample == "did"){
    .color = sym("Party")
    scale_color = c( "#0015BC", "#DE0100", "black")
    
    out0 <- .data %>% 
      filter(Party != "Difference in Party" & Wave != "During")
    
    
    # Set breaks for facts and seriousness facets
    if (str_detect(deparse(substitute(.data)), "fact")){
      .xintercept = 50
      .x = "Estimated votes (per 100 legitimate votes cast)"
      
      position_scales <- list(scale_x_continuous(breaks = seq(10, 30, 10), limits = c(5,35)),
                              scale_x_continuous(breaks = seq(10, 30, 10), limits = c(5,35)),
                              scale_x_continuous(breaks = seq(-10, 15, 5), limits = c(-10, 15)),
                            scale_x_continuous(breaks = seq(-10, 15, 5), limits = c(-10,15)))
    } else if (str_detect(deparse(substitute(.data)), "abstract")){
      .xintercept = .5
      .x = "Estimated seriousness (scaled emotional reaction score)"
      position_scales <- list(scale_x_continuous(breaks = seq(0.7, 0.9, 0.1), limits = c(0.65,0.95)),
                              scale_x_continuous(breaks = seq(0.7, 0.9, 0.1), limits = c(0.65,0.95)),
                              scale_x_continuous(breaks = seq(0, 0.1, 0.05), limits = c(-.025, 0.1)),
                              scale_x_continuous(breaks = seq(0, 0.1, 0.05), limits = c(-.025,0.1)))
    } 
    
  # Panel
  } else if (.sample == "panel"){
    .color = sym("Party")
    scale_color = c( "#0015BC", "#DE0100")
    
    out0 <- .data %>% 
      filter(term == "(Intercept)")
    
    
    # Set breaks for facts and seriousness facets
    if (str_detect(deparse(substitute(.data)), "fact")){
      .xintercept = 50
      .x = "Estimated votes (per 100 legitimate votes cast)"
      
      position_scales <- list(scale_x_continuous(breaks = seq(10, 30, 10), limits = c(5,35)),
                              scale_x_continuous(breaks = seq(10, 30, 10), limits = c(5,35)),
                              scale_x_continuous(breaks = seq(10, 30, 10), limits = c(5,35)),
                              scale_x_continuous(breaks = seq(10, 30, 10), limits = c(5,35)))
    } else if (str_detect(deparse(substitute(.data)), "abstract")){
      .xintercept = .5
      .x = "Estimated seriousness (scaled emotional reaction score)"
      position_scales <- list(scale_x_continuous(breaks = seq(0.7, 0.9, 0.1), limits = c(0.65,0.95)),
                              scale_x_continuous(breaks = seq(0.7, 0.9, 0.1), limits = c(0.65,0.95)),
                              scale_x_continuous(breaks = seq(0.7, 0.9, 0.1), limits = c(0.65,0.95)),
                              scale_x_continuous(breaks = seq(0.7, 0.9, 0.1), limits = c(0.65,0.95)))
    } else if (str_detect(deparse(substitute(.data)), "right")){
      .xintercept = 3.5
      .x = "Estimated right (scaled emotional reaction score)"
      position_scales <- list(scale_x_continuous(breaks = seq(5, 7, 0.5), limits = c(5,7)),
                              scale_x_continuous(breaks = seq(5, 7, 0.5), limits = c(5,7)),
                              scale_x_continuous(breaks = seq(5, 7, 0.5), limits = c(5,7)),
                              scale_x_continuous(breaks = seq(5, 7, 0.5), limits = c(5,7)))
    }
    
  } 
  
  # 2. Base plotting device for all samples
  g <- ggplot(out0)+ 
    geom_pointrange(aes(y = Error,
                        x = estimate,
                        xmin = conf.low,
                        xmax = conf.high,
                        color = !!.color, 
                        shape = !!.color),
                    position = position_dodge2(width = .75, reverse = T)) +
    theme_few() +
    scale_y_discrete(limits = rev) +
    labs(y = NULL,
         x = .x)  +
    guides(color = guide_legend()) +
    theme(legend.position = "bottom")
  
  # 3. Apply additional plotting depending on sample
  ## All
  if (.sample == "all_respondents"){
    g <- g + 
      geom_text(aes(y = Error,
                    x = estimate, 
                    label = format(round(estimate, 2), nsmall = 2)),
                nudge_y = .2) + 
      geom_vline(data = .data %>% filter(outcome == "abstract_values"),
                 aes(xintercept = .5),
                 linetype = "dashed", 
                 color = "gray") +
      geom_vline(data = .data %>% filter(outcome == "fact_values"),
                 aes(xintercept = 50),
                 linetype = "dashed", 
                 color = "gray")+
      theme(axis.text.y=element_text(vjust=0.5)) + 
      facet_grid(~ factor(outcome, levels = c("fact_values", "abstract_values"),
                                 labels = c("Frequency of Error", "Seriousness of Error")),
                        scales = "free_x") +
      ggh4x::facetted_pos_scales(x = position_scales)
    ## By Party
  } else if (.sample == "by_party"){
    g <- g + 
      geom_text(data = .data |> filter(Party == "Democrat"),
                aes(y = Error,
                    x = estimate, 
                    label = format(round(estimate, 2), nsmall = 2)),
                nudge_y = .35) + 
      geom_text(data = .data |> filter(Party == "Republican"),
                aes(y = Error,
                    x = estimate, 
                    label = format(round(estimate, 2), nsmall = 2)),
                nudge_y = -.35) + 
      scale_color_manual(values = scale_color) +
      
      theme(axis.text.y=element_text(vjust=0.5)) + 
      facet_grid(~ factor(outcome, levels = c("fact_values", "abstract_values"),
                          labels = c("Frequency of Error", "Seriousness of Error")),
                 scales = "free_x") +
      ggh4x::facetted_pos_scales(x = position_scales)
    
  ## For DiD figures
  } else if (.sample == "did"){
    g <- g + 
      scale_color_manual(values = scale_color) +
      geom_vline(data = .data %>% filter(Wave  %in% c("Difference in Wave", "Difference in Differences")),
                 aes(xintercept = 0),
                 linetype = "dashed",
                 color = "gray") +
      # geom_vline(data = .data %>% filter(!(Wave %in% c("Difference in Wave", "Difference in Differences"))),
      #            aes(xintercept = .xintercept),
      #            linetype = "dashed",
      #            color = "gray")+
      facet_grid(~ factor(Wave, levels = c("Pre", "Post",
                                           "Difference in Wave", "Difference in Differences"),
                          labels = c("A) Pre, 2020\n(October 2020)", 
                                     "B) Post, 2020\n(January 2021)", 
                                     "C) Post - Pre", 
                                     "D) Difference in Differences\n(Post - Pre)\n(Republican - Democrat)")),
                 scales = "free_x") +
      ggh4x::facetted_pos_scales(x = position_scales)
  
  ## Panel
  } else if (.sample == "panel"){
    g <- g + 
      scale_color_manual(values = scale_color) +
      # geom_vline(data = .data %>% filter(Wave  %in% c("Difference in Wave", "Difference in Differences")),
      #            aes(xintercept = 0),
      #            linetype = "dashed",
      #            color = "gray") +
      # geom_vline(data = .data %>% filter(!(Wave %in% c("Difference in Wave", "Difference in Differences"))),
      #            aes(xintercept = .xintercept),
      #            linetype = "dashed",
      #            color = "gray")+
      facet_grid(~ factor(Wave, levels = c("Pre", "During", "Post", "2024"),
                          labels = c("A) Pre, 2020\n(October 2020)", 
                                     "B) During, 2020\n(November 2020)",
                                     "C) Post, 2020\n(January 2021)", 
                                     "D) Pre, 2024\n(May 2024)")),
                 scales = "free_x") +
      ggh4x::facetted_pos_scales(x = position_scales)
    
  }
  
  return(g)
  
}


